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ABSTRACT 


The transport of chemically reactive solutes through porous media is 
generally affected by multiple reactions with different rate constants. In this study, 
transport phenomena affected by rate limited sorption, first-order mass transfer and first- 
order transformation is studied. The governing equations are taken in their most general 
form for the analysis. Analytical solutions are developed using the laplace transform 
approach. 


Exact solutions for batch reactor problems are arrived at. Semi-analjdical 
solutions are developed for lab scale transport systems with constant hydrodynamic 
dispersion coefficient. Analytical expressions are developed for temporal moments of the 
plume. An effective numerical technique to invert the laplace transform is used to obtain 
the break-through curves. The temporal behavior of the phenomenon under different rate- 
constant values is studied. The results firom the present study are compared with the field 
experiments and they are found matching closely. 

These solutions will help in determine a criterion for the permissible time 
step length for a good numerical solution of the transport problem and to test the 
accuracy of the simulation codes. They will also be useful as a prehminary analysis tools 
for ascertaining the relative importance of various processes for given conditions. 


VI 


TABLE OF CONTENTS 


Title 

Dedication 

Acknowledgements 

Certificate 

Abstract 

Table of contents 
Notations 
Index to tables 
Index to figures 

INTRODUCTION 

1 . 1 General Introduction 

1 .2 Review of Literature 

1 .3 Objective of the present study 

MODEL 

2.1 Introduction 

2.2 Model Description 

2.3 Governing Equations 

MATHEMATICAL ANALYSIS 


3.1 Introduction 



3 .2 Assumptions made about the solute and medium 

3.3 Nondimensionalization 

3.4 Transforming into laplace domain 


10 

11 

14 


SOLUTIONS 


4.1 

Introduction 

16 

4.2 

Batch Reactor System 

16 

4.3 

Column Transport Systems 

21 


4.3.1 Analytical Solution 

21 


4.3.2 Numerical Inversion 

25 


4.3.2. 1 Introduction to the method 

25 


4.3. 2.2 Methodology 

25 


4.3. 2.3 Application to the present problem 

28 

RESULTS AND DISCUSSIONS 


5.1 

Introduction 

29 

5.2 

Data used 

29 

5.3 

Batch Reactor System 

30 

5.4 

Column Systems 

31 


REFERENCES 43 



NOTATIONS 


Ca 

c„ 

Co 

CaO 

CnO 

D 

dW 

f 

F 

G(p) 

H(p) 

H(0) 

K 

ka2 

kn2 

1 

m 

M„‘ 

P 

q 

Sa2 


Input parameter in inversion technique 
Concentration of solute in advective region 
Concentration of solute in nonadvective region 
Reference concentration 

Initial concentration of solid phase solute in advective region 

Initial concentration of solid phase solute in nonadvective region 

Hydrodynamic dispersion coefficient 

Effective dispersion coefficient 

Mass fraction of the sorbent in the advective region 

Fraction of sorbent with instantaneous sorption 

The function of the laplace variable with a ratio of polynomials 

1/G(p) 

H(p)lp=o 

Equilibrium sorption coefficient 

First order reverse sorption rate-constant in advective region 

First order reverse sorption rate-constant in nonadvective region 

Input parameter in the inversion technique 

Input parameter in the inversion technique 

n* absolute temporal moment of the concentration 

Laplace variable 

Specific discharge 

Rate-limited sorbed phase concentration in advective region 


IX 



SaO 

Sn2 

SnO 

t 


to 

T 


X 

Ga 


P 

a 


Pal 

Pa2 


Initial concentration of the sorbed phase solute in advective region 
Rate-limited sorbed phase concentration in nonadvective region 
Initial concentration of sorbed phase solute in nonadvective region 
Time 

Pulse duration 
Reference time 
Effective velocity 
Distance 

Fractional volumetric water content of advective region 
Fractional volumetric water content of nonadvective region 
Bulk density of the porous medium 
First order mass transfer coefficient 
Transformation coefficient for instantaneous sorption 
Transformation coefficient for rate-limited sorption 
n*^ absolute moment of Ca normalized with mass 


X 



INDEX TO TABLES 


3. 1 Nondimensional parameters for advective region 1 1 

3.2 Nondimensional parameters for nonadvective region 12 

4. 1 Forms of G(p) for various levels of complexity in models 21 

4.2 Contributions of various terms in Ca towards its time domain equivalent 24 

5. 1 Data used for the present study 30 

5.2 Errors with various parameter values for inversion 34 

5.3 Temporal moments for various values of k 34 

5.4 T emporal moments for various values of a 35 


XI 



INDEX TO FIGURES 


5 . 1 Effect of k on concentration 3 6 

5.2 Effect of a on concentration 36 

5.3 Concentration profile at quarter length of the column 37 

5 .4 Concentration profile at half length of the column 3 7 

5.5 Concentration profile at three-fourth length of the column 38 

5.6 The break-through curve 38 

5.7 Effect ofk on break-through curve 39 

5.8 Effect of a on break-through curve 39 

5.9 Effect ofk on the skewness of the break-through curve 40 

5.10 Effect of a on the skewness of the break-through curve 40 

5.11 Effect of k on the second moment 4 1 

I 

5.12 Effect of a on the second moment 41 

5.13 Concentration profile at mean arrival time 42 

5.14 Comparison with observed data 42 


XU 



Chapter 1 

INTRODUCTION 


1.1 General Introduction 


The quality of water is of major interest in any development and 
management of water resources system. With the increased demand for water in most 
parts of the world, and with the intensification of water utilization, the quality problem 
becomes the limiting factor m the development of water resources. The quality of both 
surface and groundwater resources is deterioring as a result of pollution. Special attention 
should be devoted to the pollution of ground water due to their very slow velocities. 
Although it' seems that groundwater is more protected than surface watK' against 
pollution, it is still subject to pollution, and when the contamination of groundwater 
occurs, the restoration to the original, non-polluted state is more difficult. Due to the 
potential of the problem, the accuracy in the estimation of the extent of pollution assumes 
great importance. In technical terms, this problem of groundwater pollution is treated as 
mass transport phenomenon in porous media, where the considered mass is that of some 
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solute or pollutant moving with the solvent which is water in the interstices of the porous 
medium. 


The principle mechanisms affecting the transport of a solute in a porous 
medium are convection, mechanical dispersion, molecular diffusion, solid-solute 
interactions and various chemical reactions and decay phenomena, which may be 
regarded as source-sink phenomena for the solute. 

For mass transport phenomena, mathematical analysis and quantification 
starts at mass conservation equation which states that the difference between mass inflow 
rate and mass outflow rate should be equal to the change in mass storage with time. For 
systems, which involve processes like advection, diffusion, dispersion and etc, and" if the 
tracer is ideal, this expression suffices. 

But when the phenomenon involves a chemically reactive element in the 
mass that is being transported, the mass conservation equation should be expanded so as 
to incorporate this also. The material, while getting transported, may disintegrate or the 
porous medium, depending on its composition, may contribute to the mass. The mass 
conservation equations gets modified in terms of source or sink terms, depending on 
whether a constituent is being added to the mass or removed by chemical processes. Then 
the appropriate statement for mass conservation when reactions are considered would be. 

Mass inflow rate - Mass outflow rate ± Mass production rate = Change of mass with time. 

The plus or minus sign designates either a source or a si^ and takes on 
different forms for different reactions. For convenience, the somce or sink terms are 
quantified as mass, produced or consumed per unit volume per unit time. Transport 
involving several dissolved species in the transport requires a system of mass conservation 
equations, one for each constituent. Appropriate rate laws are used to express the reactions 
mathematically, depending on the reaction. 
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Reactions are classified as homogeneous, operating within a single phase, 
or heterogeneous, operating between two phases, such as, solid and liquid. Also, reactions 
can be described from an equilibrium or kinetic point of view depending on the rate of the 
reaction relative to the mass transport process. 

1.2 Review of Literature 


The transport of chemically and biologically reactive solute through 
porous media, is affected by various rate-limited processes (Brusseau, 1989). Several 
techniques have been reported in the literature for numerical simulation of the advective- 
dispersive-reactive (ADR) equation in the presence of rate-limited processes. .These 
include the modified method of characteristics and mixed finite elements method (Chiang 
et al 1989), split-operator, Petrov-Galerkin methods (Miller, 1993). Fry [1993] reported 
an analytical solution to the solute transport equation with rate-limited desorption and 
decay. Tompson [1993] gave a numerical simulation technique for chemical migration in 
physically and chemically heterogeneous porous media. Daus et al [1985], tested various 
finite element formulations of the advection-dispersion equation and gave a comparative 
error analysis among those models. Relative errors are tabulated for various different 
models. An analytical solution for describing the transport phenomena m heterogeneous 
media with a distance dependent dispersion relation has been developed by Yates [1990]. 
He, first transformed the advection-dispersion equation with space-dependent dispersion 
coefficient, in to laplace domain, and then used contour integration techniques to invert in 
to original time domain. Goltz et al [1986], solved the three dimensional solute transport 
problem in infinite medium with mobile and immobile zones. Brusseau et al [1992], 
discussed the difficulties and factors to be considered while modeling transport 
phenomenon influenced by multi-process nonequilibruim and transformation reactions. 
Valochchi [1992], discussed the accuracy of operator splitting techniques of ADR 
problems. 
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The formulation suggested by Brusseau [1991] is quite general and 
includes various other models as its special cases. It is therefore used as the basis for the 
present study. 


1.3 Objective of the present study 


With the advent of high speed computing facilities, numerical methods are 
gaining more popularity for modeling reactive chemical migration through porous media. 
In such cases, evaluating the effect of multiple reactions on the time discretization 
becomes critical. The objective of the present study is to provide a basis to address that 
problem, by developing analytical and semi-analytical solutions for one-dimensional 
transport of reactive contaminants through porous media undergoing rate-limited sorption 
and first order mass transfer and transformation at a generic level. Laplace transform is 
used to solve the equations. The batch reactor systems or point-process systems and 
laboratory scale column transport systems are to be solved analytically, at a generic level. 
Expressions are to be developed to obtain the quantities of practical significance and 
firom them, the temporal dependency of the phenomenon is to be studied. The analytical 
solutions are to be tested by using numerical techniques to invert the equations in laplace 
domain in to original time domain. The effect of the rate constants on the transport 
phenomenon is to be studied. The results from the study are to be compared with the 
experimental data available in the literature. 
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Chapter 2 
MODEL 


2.1 Introduction 

Theories on dispersion, so far published, suggest two approaches to model 
mass transport phenomenon in porous medium. The main effort has been to express 
hydrodynamic dispersion macroscopically through a partial differential equation and to 
determine the nature of coefficients, which appear in this equation. Literature suggests 
two broad ways to model this problem. In one approach, the porous medium is replaced 
by a fictitious, greatly simplified, model in which the spreading of a solute can be 
analyzed by exact mathematical methods. Single capillary tubes, a bxmdle of capillaries, 
an array of mixing cells, are examples of such models. The second approach is to 
construct a conceptual statistical model of the microscopic motion of solute particles and 
to average these motions in qrder to obtain a macroscopic description of them. 

Many types of models used to describe the transport of solutes in porous 
media are based on the convection-dispersion equation as given by Bear [1972, 1979]. 
The convection-dispersion equation may include terms, which describe the hydrodynamic 
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dispersion, fluid advection, linear equilibrium, first-order reaction and possibly zeroth 
order production. Models describing the chemical reactivity of the contaminant and the 
phase, in which it occurs, also use the same model. 

The reactions can be categorized as biodegradation, radioactive decay, 
hydrolysis etc. In general, models make assumptions about the phase of the contaminant, 
whether the contaminant exists in solid phase or hquid phase or both. The mixing is 
assumed to occur in various domains, such as advective, nonadvective. The complexity 
of the model is influenced by the occurrence of mass transformation and rate limited 
sorption in various phases. Allowing these reactions makes the model more complex. 

2.2 Model Description 

For the present study, a two-phase, two-domain contaminant transport 
model has been assumed. Apart from mechanical mixing and microdispersion, the 
transport is influenced by reactivity of the solute. In this model, the entire flow domain is 
assumed to be divided into two regions namely, 1) Advective Region 2) Nonadvective 
region. The mass balance equation is applied on both regions separately. At the grain 
level, the mixing is influenced by the sorption on the surface of the particle. Every 
particle in the flow domain consists of two sites on surface namely equilibrium site and 
rate limited site. The contaminant is assumed to be chemically reactive and to be existing 
in both solution phase and solid phase. The rate-limited sorption is used to represent 
chemical reaction and diffusional processes and the first-order mass transfer is used to 
represent the macro-scale heterogeneity by mass transfer between advective and non- 
advective regions of any size, shape or type. The transformation coefficients are used to 
represent irreversible transformation reactions at a generic level. First order 
transformation accurately describes the radioactive decay. Biodegradation is usually 
specified to occur in the solution phase but radioactive decay and hydrolysis reactions 
may occur in any phase. That means, transformation is allowed to occur in any phase and 
any domain. 
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2.3 Governing Equations 


The governing equation for the transport of a reactive chemical through 
the advective region of a porous medium is given by (Brusseau, 1991) 




bd: 


. ac* 1 . dc: 


dx 




j j 


dx] (2.1) 


- « ’ ( c: - c; ) - ( e. + / p )c: - m;. / p 


Where, 

6a is fractional volumetric water content of the advective region (L*’), 
f is mass fraction of the sorbent constituting the advective region (M°), 
p is the bulk density of the porous medium (ML‘^), 

F is the fraction of the sorbent with instantaneous sorption, 

K is the equilibrium sorption coefficient (L^M"'), 

C is the concentration of the solute in solution (ML'^), 
t is the time (T), 

S 2 is the rate- limi ted sorbed-phase concentration (M hT^), 

Xi is the coordinate directions (L), 

Dij is the hydrodynamic dispersion tensor (L^T’^), 
qi are the specific discharge components (LT^), 

a is the first order coefficient for^mass transfer between the advective and 
nonadvective regions (T*), 

Pi and P 2 are the transformation coefficients (T'^) for the instantaneous and 
rate-limited sorbed-phase domains, respectively. 

Subscript a, refers to the advective region. Similar parameters in the 
nonadvective domain are denoted with the subscript n, and the superscript * is used to 
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distinguish the notations for these variables from their nondimensional counterparts 
which will come in the later sections. 

The mass balance of the nonadvective domain can be described by the 
following equation. 


(e, + (i-/)pF.ji:.)^+(i-/)p^ = a[c:-c:) 

- ( 0, + m: (1 - /) p f, k, )c; + (1 - /) p s; 


( 2 . 2 ) 


Dynamics of sorption and transformation for the rate-limited domains are 
described by the following equations. 


a,- ‘ 


- 

(2.3) 

= k- [( 

dt‘ " 


- 

(2.4) 


Where, k*a2 and k*n2 are the first-order reverse sorption rate coefficients 
(T’^) for the advective and nonadvective regions, respectively. 


The laplace transform can be applied on a fimction so as to make the 
behavior of the function smoother than in the original domain. 'LG\f(t) is a real valued 
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function defined over the interval (-a, a) such that f(t) - 0 for all t < 0. The laplace 
transform of f(t), denoted by 'Lifft)} is defined as, 

f(t)e-ds (2.5) 

Where, parameter s is called laplace variable, independent of space and 
time and can be real or complex number. The laplace transform of f(t) is said to exist if 
the integral (2.5) is convergent for some value of s. 
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Chapter 3 

MATHEMATICAL ANALYSIS 

3.1 Introduction 

\ 

This chapter describes the mathematical analysis done to solve ADR 
equation analytically. The governing equations are expressed according to the 
assumptions made regarding the solute and medium. Then the equations are 
nondimensionalized using appropriate nondimensional parameters. The equations in the 
nondimensional form are transformed in to the laplace domain using the definition of 
laplace transform, given in section 2.3. 

i 

3.2 Assumptions made about the .solute, and medium 

In order to solve the one-dimensional classic ADR equation anal 3 dically, 
certain assumptions are made regarding the phenomenon of solute transport in porous 
medium. They are, 

1) The medium is homogenous and isotropic. 

2) The phenomenon is one-dimensional. 
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3) The medium is fully saturated. 

4) Transformation coefficients are used to represent irreversible reactions. 

3.3 Non-Dimensionalization 


The governing equations are non-dimensionalized so as to get the 
equations in a form free from the units of the parameters being used for calculations. The 
equations are made more tractable to solve anal 54 ically, by doing so. The non- 
dimensional parameters are defined for advective domain and non-advective domain 
similarly. For advective domain, they are given in table 3.1. 


Table 3.1 Nondimensional parameters for advective region 


Parameter 

Definition 

Ca 

CVCo 

Sa 

SV(l-Fa) KaCo 

t 

tVT 

D 

BaD'e Ra,/q’^T 

X 

x*0 Rat/q*T 

Ral 

(t)a+fp(l-Fa)Ka/e 

Ra2 

fpFaKa/e 

(0 

a*T/9 

4al 

( ll*a't>a + p‘aifpFaKa/9 ) T 


Ra2 l-t a^T 

ka2 

k*a2T 

Va 

ka2 31 ^ 

Yai 

RaiVa + Raskai 

^ a2 : 

Va^l "^ka2^al 


SWTRAL LiBRAK' 

I. i.T., KANrUt 
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Where, G = 0a + ©„ 

<})a = 0a/0 


Sai=KaCa 

Sn.=K„C„ 


Sal 3nd Sni are the concentrations of sohd phase solute in the equilibrium 
sites of advective and nonadvective regions, respectively. Similar nondimensional 
parameters are defined for the nonadvective domain variables also. They are given in 
table 3.2. 


Table 3.2 Nondimensional parameters for nonadvective region 


Parameter 

Definition 

Cn 

CVCo 

Sn 

S*m/(l-F„)K«Co 

Rnl 

(t)„+(l-f)p(l-F„)Kn/e 


(l-f)p(l-F„)K„/e 

4nl 

( P*n<|)n + p*ni(l-f)pF„Kn/0 ) T 


Rci 2 fJ' n2T 

km 

kmT 

Vn 

km+ P*n2T 

Ynl 

RfilVn ^nl Rn2kii2 

2 

Y 

Vn^nl kn2^nl 


Where, <})„ = 6n/0 

The above-defined nondimensional parameters are substituted in the 
equations, (2.1), (2.2), (2.3) and (2.4) and are simplified. It is assumed that the volumetric 
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water content and the dispersion coefficient are spatially invariable. Also for notational 
convenience, it is assumed that neither Rai nor Rm are zero. If Rni is zero, all the 
parameters in the nonadvective domam can be set to zero. The mass balance equations 
for advective and nonadvective domain takes the form as given below, respectively. 


R I R = R D^-^-R (r -C r 

(3.1) 


i?„i — + R 
dt 


nl 


dS 

n 

dt 


= 0 ) 




(3.2) 


The governing equations for the dynanaics of rate limited sorption for both 
advective and nonadvective domain are also nondimensionalized. They take the form 
given below. For the rate-limited sorption in advective domain, 


5 ^ 

dt 


= K2C,-^Ja 


(3.3) 


The governing equation for rate-limited sorption in nonadvective domain, 
after nondimensionalization, 


dt 




(3.4) 


The equations, (3.1) and (3.2) together with the sorption dynamics 
equations, (3.3) and (3.4) can be solved to obtain the non-dimensional concentrations Ca, 
Cn, Sa 2 and Sn 2 for given initial and boundary conditions. 
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3.4 Transforming into laplace domain 


The equations given in the previous section involve partial derivatives and 
solving it analytically, in real domain needs many assumptions. To solve them in their 
existing form, they are transformed into laplace domain using the definition given in 
(2.5). The laplace transform is used to eliminate time dependency in the equations. In 
laplace domain, solutions are generally smoother than in the time domain. In addition, the 
solutions for different times are independent (or parallel), in the sense that, solution for 
one time is not dependent on solution for some other time, which means that the problem 
of numerical error accumulation is not encountered. 


Taking the laplace transform of the equations (3.2)-(3.4), and assuming 
that the initial mass of solute can exist in any phase and any domain, yields the following 
expressions. The rate-limited sorption equations become. 


^ _ ^a2^a '^^aO 


(3.5) 


P+^n 


(3.6) 


The mass balance equation for nonadvective domain becomes. 


C„ = 


— + 


R 


nl 




’nO 


Yr,+(0 


r„+® 


Yn+G) 


^nl ^nl 




(3.7) 


In which, 


Yn 


RnlP^+YnxP^yl2 

P+^n 


(3.8) 
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Where, p is the lapalace variable, which is independent of space and time 
and can be real or complex. The overbar indicates the concentration equivalents in 
laplace domain. Sao is the initial concentration of solute in sohd phase and in advective 
domain, Sno, in solid phase and in nonadvective domain, and Cao is the concentration of 
solute in advective domain, liquid phase. The subscript a, denotes advective domain and 
n, nonadvective domain. In order to get the laplace equivalent of Ca, the equation (3.1) 
would be transformed in to laplace domain depending on the assumed conditions about 
the phenomena. They will be described in the next chapter. 
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Chapter 4 
SOLUTIONS 


4.1 Introduction 

This chapter describes the solutions for two different transport systems 
namely batch reactor system and column transport system. Laplace transform is applied 
on equation (3.1) for both the cases separately, depending on the assumptions made. Then 
they are inverted from laplace domain to original time domain, in order to obtain the 
hondimensional concentrations. 

4.2 Batch Reactor System 


For batch reactor problems, the phenomenon is assumed to be local. That 
means the process is essentially treated as point process. The concentration is considered 
to be spatially invariable. With this assumption, the spatial terms in the equation (3.1) 
will become zero. Also, for point processes, the reference time T is taken as one unit. 
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Applying laplace transform on (3.1) with the above assumptions and 
substituting (3.2), (3.2) and (3.4) results in the following expression. 


7^+0)- 


y« 


y„+a) 


-Kl^aO + 


y„+co 


C 


nO 


I ^a2^a2 n , 


P+^a 


CO 

^n2 ^n2 

_r„+®_ 

.P-^^n _ 


S„o (4-1) 


Where, Ya is given by. 


^ K^p"+ra,p+rl2 

* ^ 

Denoting the term in the square bracket on the left hand side by G^), 
equation (4. 1) is represented as. 


1 

G(/>) 



R., C.„ + 


7„+0) p+v^ 


aO 


+ 


CO 

^n2^n2 

_7„+CD_ 

.P-^^n _ 


. (4.2) 


For two domain models, R„i is not equal to zero. The function G(p) in its 
most general form is expressed as a ratio of two polynomials as given. 


G{p) = 


hp" +Kp^ 

a,p* +a^p^ +a,p^ +a,p+a. 


(4.3) 


For single domain models, Rni is zero. Then the function Gip) is given by. 


G(p) = — 
r» 


P + ^a 

Ra\P^+ralP+Yl2 
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Where, 


bl =rll+yn\^a^G)^n 
bQ=yl2V^+(OV^V^ 


a, =r.,^„.+r„l^al+®^al +<»^nl 

=^«ir;2 +Kxyl2 +®r„, +®y.: +r.i7„i 

^1 = r„l 7^2 + 7«1 7'2 + <»7a2 + ®7'2 + mal^n + 

«o = 7^2 yli + «07 '2 + (oyl2 


Now, G(p) is taken to the right hand side and the coefficients of Cao, Qio, 
Sao and Sno are multiplied with the same. This results in an expression of the form, 


c. =G„(/>)C„ +G^(p)C., +G„(p)S.. +G„(p)S., (4.4) 


Where, all the G’s are functions of p and are expressed similar to that of 
G(p). Denominators of each are same as that of G(p). The expressions of numerator 
differ. 


For the most general case, coefficients in the numerator of Gac(p) are 

given by. 
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^3 =Rar^n^ 

^1 = ^alY^ +Kjnl^a 

h =Kjn2^a+K,(0^a^n 


The Coefficients in the numerator of G„c^) are given by. 


b,=R„,co 

b^=R„^cov^+R„,cov^ 


The coefficients in the numerator of Gas(p) are given by, 


b 2 =RaxRnxK 2 

^:=^<,2^.2r„:+^„2^.2® 


Similarly, G„s(p) is given by. 


i>X =^„2^n2® 

K=Rn2K2<^^a 
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After knowing all the coefficients, the equation (4.4) can be solved using 
simple mathematical methods. The roots of the polynomial in the denominator can be 
found by standard anal 3 dical methods available for solving a biquadratic equation 
(Abramowitz and Stegun, 1970). After finding all the zeroes in the denominator, GQj) can 
be expressed as partial fractions. 


G(p) = 


i=i p+p. 


(4.5) 


Where, pi are all real and noimegative and are the negatives of the zeroes 
of the denominator, P/ are the coefficients of the partial fractions, and N is the total 
number of roots, which may vary from one to four depending on the type of the model. 

Substituting all the G(p)s in the equation (4.4) in their partial fraction form 
and inverting in to the original time domain gives a relation, as given below. 


II 


o 

+ 







L.-=i j 




_i=l 




(4.6) 


The above expression describes the behavior of nondimensional 
concentration in time domain. The temporal variations can be obtained by substituting the 
value of time t. The nondimensional concentrations in other regions and phases can be 
obtained using the equations, (3.5), (3.6) and (3.7). Though, the e^qjression given above 
for G(p), is the most general one, it changes depending on the type of the model. For various 
complexity levels, G(p) has been arrived at. They are givai in Table (4.1). 
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Table 4.3 Forms of G(p) for various complexity levels in models 


CASE 

MODEL 

G(P) 

1 

Single Region, LEA 

1 /(p+a3) 

2 

Single Region, Kinetics 

(p + bi) / (p^ + asp + a2) 

3 

Double Region, LEA 

(p + b2) / (p^ + asp + as) 

4 

Double Region, Kinetics 

(b3p^+b2p^+bip+bo)/(a4p'‘+a3p^+a2p^+aip+ao) 


4.3 Column Transport Systems 


4.3.1 Analytical Solution 


For column transport systems of laboratory scale, and in experiments with 
packed columns, the porous medium can be assumed to be homogeneous. The 
dispersivity can be taken as constant in space and time. Though, with this assumption, 
most of the field-scale heterogeneity is not represented, estimations with these 
simplifications give reasonable results. The reference length L, is taken equal to the 
length of the column and the reference time T is. 



Using the above assumptions, equation (3.1) is transformed into laplace 
domain. The resulting expression is, 


1 

R.Mp) 



SC^ 

dx 


(4.7) 
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In order to solve the above equation, a constant boundary flux condition 
has been used. Initially, there will be no solute in the column and a pulse of unit 
concentration, will be introduced for certain time and after that time, incoming 
concentration will be zero. That means, if to is the dirration of the pulse, the boundary 
condition at the upstream would be. 


C=lforO<t<to 
C = 0 for t > to 

This type of boundary condition is called first type boxmdary condition. 
The downstream boundary condition is, concentration gradient becomes zero at a large 
distance. With these conditions, solution of these type of equations is given by Van 
Genuchten [1980], and is. 


C = 




■e- 


1-J1+ 


AD 


RalGip) 


X 

2D 


( 4 . 8 ) 


Analytical inversion of this equation to the time domain is not possible for 
general cases. It can be inverted analytically for very simple cases of transport only, for 
example, transport with no rate-limited sorption, dual porosity or degradation. In such 
cases, G(p) reduces to, G(p) = 1/p. Therefore, a numerical inversion method has been 
used for the pmpose and would be explained later in this section. 

For column transport systems, break-through curves are generally more 
useful than the spatial distributions. The temporal moments of the plume at any distance 
can be obtained using the Aris method, as used by Goltz et. al., [1987] and Valocchi 
[1990]. In this method, the exponential term in the definition of laplace transform is 
expanded m the powers of (pt) and are differentiated successively. This results in. 
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(4.9) 



/? = 0 


Where, Mn* is the nth absolute temporal moment of Ca. 

The quantity of practical significance is the first absolute moment, 
normalized with mass, which is nothing but the mean arrival time of the break-through 
curve. Also, the normalized second and third moments taken about the mean arrival time, 
are needed to estimate physical parameters like dispersion and skewness. 

Using (4.9) directly will make the evaluation of the same difficult, as the 
expression (4.8) involves exponential terms, and also is a ratio. Instead, a slightly 
modified version of (4.9) will be used, which will directly give normalized first moment 
and the normalized central second and third moments. It is given as, 




dp^ 


p=0 


(4.10) 


Where, pn is the nth absolute temporal moment of Ca normalized with 
mass. The above equation has an advantage that any factor involving exponential is 
greatly simplified. Also, the transformed solution would involve multiplication of several 
factors, which can be evaluated separately. 

Using this approach, the contributions of various terms in the equation 
(4.10) towards the mean arrival time, the variance, and the skewness are listed in the table 
(4.2), given below. 
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Table 4.4 Contributions of various terms in Ca towards its time domain equivalent 


Fimction 

(l-exp(-pto))/p 

Exp{-0.5x [V(l+4DH(p)/Rai)-l]/D} 

Mo* 

to 

Exp(-x(Hi-l)/2D) 

/ 


to/2 

XH’(0)/Hi 

^2* 

to^/12 

2DxH’(0)VHy - xH”(0)/Hi 

143* 

0 

12D^xH’(0)VHy - 6DxH’(0)H”(0)/Hy + xH”’(0)/Hi 


Where, 

H(p)=l/G(^; 
Hi=V(l+4DH(0)) 
H2 = Hi(1 + Hi) 


The expressions for H(0) and its successive derivatives are given as, 

H(0) = ao/bo 


H’(0) = (aibo-aobi)V 


H”(0) = 2 (aobi^ - aobob2 + a2bo^ - aibobi)/bo^ 


(4.11) 


H”’(0) = 6 (2aobobib2 - aibo^b2 - aobi^ + aibobi^ - a2bo^bi + asbo^ - aobo^ybo"* 

The above-derived moments can be used to find out the physical quantities 
like velocity and dispersion coefficient. The effective velocity of the flow is given by. 



(4.12) 


24 















According to the fickian model of dispersion, dispersion coefficient can be 
obtained from the break-through curves. The effective dispersion coefficient can be 
related to the moments of the break-through curve using the expression, 




( 4 . 13 ) 


4.3.2 Numerical Inversion 


4.3.2.1 Introduction to the method 


In order to invert the equation (4.8), which is in the laplace domain, in to 
the original time domain, a numerical technique has been used. Hosono [1981], suggested 
a method for ILT problems. The essential point of this method is to approximate the 
exponential term in the definition of ILT by the function, Eac = exp(a)/cosh(a-s). This can 
be expanded in to an infini te series from which an approximate of inverse transform can 
be obtained. 

4.3.2.2 Methodology 

The inverse laplace transform is defined as, 

1 y+y® 

r ‘ {F(i)} = f(t) = -- KC-s) e” * ('t- 14) 

2>Fr-> 
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Let, 1) F(s) is defined in Re s > 0 


2) F(s) is regular 

3) I s 1 — ► 00 , F(s) — ► 0 

afc ]|e ^ 

4) F (s) = F(s ), where, is the complex conjugate. 




eV 

/2cosh.(a-s) 

~^] f 7(-l)"/ , 

2 As-a-j(n-0.5)K] 




(4.15) 


The approximated inverse transform is defined as. 


1 r+y“ 

Ac («.<") = |F(s)£.^(i<,a)* 


On substituting (4.15) in the above expression, f.c(t, a) becomes. 


= — [^i+^2+^3+ ] (4.16) 

Where, 

F„ =(-!)” + (4.17) 
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The expression (4.16) involves an infinite series. In order to truncate the 
this infinite series, a criterion based on euler transformation has been used, which works 
under the following conditions. 

1) There exists an integer k > 1 so that the signs of Fn alternate for n^> k. 

2) Forn> k, */ 2 < | Fn+i /F„ | < 1 
Under these conditions, (4.15) becomes. 


fee (^5 


rn 


(4.18) 


Where, Amn is recursively defined by. 


A = 1 


^mn-l ^mn 


^m + 0 
{ n J 


(4.19) 


Thus, (1, m)th approximation is found by. 




/-I m 


n=l 


n=0 


(4.20) 


This expression can be used to find out an approximation for ILT 
problems. The author did not mention any criterion for deciding the input parameters 1, m 
and a. These parameters are to be decided after thoroughly testing the procedure for a 
given problem. 
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43 . 2.3 Application to the present problem 


Using the expressions (4.20), (4.19) and (4.17), a computer routine in C 
has been developed to find Ca(t) from the expression (4.8). The code, when run, prompts 
the user for the values of 1, m and a and outputs the values of time and Ca- Break-through 
curves are built using this routine and are used to find normalized first moment and 
normalized central second and third moments. They are compared with the analytical 
solutions given in the previous subsection. 


As seen in the equation (4.20), the approximate of the inverted value is a 
function of a. Hence the value of a, should be chosen carefully. The behavior of the 
function, to be inverted, in the complex domain should be understood thorou^y. As the 
method does not provide any special primitives for singularities implicitly, all the 
singularities of the function are to be avoided before applying the method. The author of 
the method mentioned that when the value of a, is taken to be 3, the relative error is 
bound by 0.25%. With the functions, whose exact inversions are available, our 
experimentation with the method gave satisfactory results with m equal to 10. A value of 
40 for 1 would suffice for the purpose. The effect of a on the inverted value and the 
relative errors are discussed in detail in the next chapter. 



Chapter 5 

RESULTS AND DISCUSSIONS 


5.1 Introduction 


This chapter presents the results obtained from the expressions developed. 
For batch reactor systems, temporal variation of solution phase concentration in the 
advective region is presented, along with a sensitivity test of the parameters, namely rate 
constant and mass transfer coefficient. For column transport systems, break-through 
curves and moments from both break-through curves and analytical expressions 
developed in the previous chapter, are presented. 

5.2 Data used 


For the present study, the physical and chemical parameters are so chosen 
so as to be representative of organic contaminant transport in a sandy aquifer (Roberts et 
al 1986; Brusseau, 1992). The data is given in table 5.1. 
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Table 5.5 Data used for the present study 


Parameter 

Value 

f 

0.9 

P 

1.81 g/m^ 

0a 

0.29 

0n 

0.04 

Fa = Fn 

0.1 

II 

0.9 cmVg 

•a 

11 

*^4 

0.3 d'* 

a 

0.04 d*' 

l-W “ l-kl “ 1-Wi2 ^ “ Mill ~ M'n2 

Od-’ 


5.3 Batch Reactor System 


The nondimensional solution phase concentration in the advective region 
of the flow domain can be obtained as a function of time using the expression (4.6). The 
concentrations in the other region and phases can be obtained using the expressions (3.5), 
(3.6) and (3.7). The most general case of two regions, nonequihbrium model has been 
used here. It is assumed that initially the solute exists only in the solution phase and in the 
advective region and is taken as umty. (i.e., Cao “ !)• Along with this assumption, the data 
given in the table 5.1 is used to arrive at the nondimensional Ca at various time values. 
The rate coefficient k*a 2 is then varied by an order of magnitude to perform sensitivity 
analysis. The resulting plot is given in figure 5.1. 

From the figure 5.1, it can be seen that a smaller value of rate constant 
leads to a slower variation of concentration with time. As the magnitude of the rate 
constant is increased, the concentration varies rapidly. A similar sensitivity test is earned 
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out on the mass transfer coefficient a also. Due to the very small value of nonadvective 
region porosity, the effect of or on the concentration profile is minimal. 


5.4 Column systems 


For column systems also, the data used is same as those used for batch 
reactor problem, which are given in table 5.1. The specific discharge q* and dispersion 
coefficient D* are chosen to be 0.03m/day and O.OOSmVday, respectively. The length of 
the column is taken as 2.0m. The column is fed with a solute pulse of unit duration, (i.e., 
to = l-O)- 


Using the expressions given in table 4.2 and expressions (4.11), the 
moments of the break-through curve are obtained. Also, the expression (4.8) is inverted 
in to the original time domain using the methodology explained in the section 4.3.2. By 
inverting, the concentration versus time plots are build from which the moments are 
found by numerical integration. The formulae used for the computation are given below. 






At 


Ml =Z 

i=l ^ 


At 


t,- 





Cai+^ai+l 


At 


1=1 


At t 
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/=i 2 


At 


At . 


-|3 



At the initial stages of arriving at the results, the moments computed using 
the above formulae are compared with their analytical counteiparts. The relative errors 
between the analytical solutions and the computed solutions are computed using the error 
definition, 


_ Analytical - Computed 

Error = x 100 % 

Analytical 


The error values so computed are used to fix the parameters required for 
the numerical inversion technique. The errors for various values of parameters 1, m and a 
are listed in the table 5.2. The value of a specific parameter is changed by keeping the 
other two parameters equal to the values suggested by the author. For subsequent trials, 
the parameter values obtained so far are used. When the value of 1 is more than 30, the 
error is very low. It indicates that for the present problem, the technique has to be used 
with a minimum 1 value of 30. By fixing the value of 1 at 40, m is varied between 10 and 
20. For any value of m greater than 10 error did not show any variation. This implies that 
m can be taken as 10. By keeping 1 equal to 40 and m equal to 10, a is changed between 2 
to 5. The error steadily decreased as the value of a is increased. For a, equal to 4, the 
percentage error in zeroth moment is zero up to the second decimal. So, for the present 
problem the values of 1, m and a are chosen to be 40, 10 and 4 respectively. For all the 
subsequent computations these values are used while using the computer routine. 

The development of concentration profiles is observed by plotting 
concentration versus time at various values of x. The time taken to reach the maximum 
concentration increased with x. The magnitude of maximum concentration decreased 
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with increasing x. This agrees with the expected behavior of tiie column system. The 
plots are given in figures (5.3), (5.4) and (5.5). The break-through curve, for which x is 
equal to the length of the column L, is shown in the figure (5.6). 

The sensitivity tests, similar to the test done for batch reactor problem for 
the parameters k and a, is performed. The corresponding break-through curves are shown 
in figures (5.7) and (5.8). The parameter imder consideration has been changed by an 
order while keeping the other parameter same as given in table (5.1). For high values of 
rate-constant k, which denotes that the reactions are close to equilibrium, break-through 
curve has very small skewness since the dispersion is very small. The maximum 
concentration is higher than that for smaller values of k. When the rate-constant is 
smaller, the break-through curve is more skewed and spanned for longer time. This infers 
that the reactions are in nonequilibruim. The corresponding temporal moments are given 
in table (5.3). When the first order mass transfer coefficient a is varied by order, the 
break-through curves does not vary much. This indicates that a has very little effect on 
the phenomenon owing to the fact that the porosity value for nonadvective region is 
small. The corresponding temporal moments of the break-through curves are given in 
table (5.4). The skewness of the break-through curves for various values of k and a are 
illustrated in the figures (5.9), (5.10), respectively. 

As noted by Goltz and Roberts [1987], in the absence of transformation, 
the velocity is independent of the rate-constant. The coefficient of skewness decreases 
with the increase in the magnitude of the rate-constant. The present study matches with 
these facts. Analytical moments matched with the computed moments with relative error 
less than 0.01%. This infers that the method for numerical inversion of laplace transform 
used in this work can be used to get a very good approximation of the ILT. 

Expressions developed in the present study are used to check with the field 
experiments. The field data illustrated by Brusseau et al [1992] is used for the purpose. 
Figure 5.14 gives the comparison. It can be seen that the results firom the present study 
matches with the field observations very closely. 
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Table 5.6 Errors with various parameter values for the inversion 


1 

m 

a 



Error in 

10 



0.176 

3.811 

87.399 

20 

6 

2 

0.112 

1.656 

7.093 

30 



0.110 

1.497 

3.225 

40 



0.110 

1.471 

3.192 


10 


0.110 

1.461 

6.786 

40 

14 

2 

0.110 

1.461 

6.786 


18 


0.110 

1.461 

6.786 



2.5 

0.116 

0.236 

1.286 



3.0 

0.087 

0.083 

0.256 



3.5 

0.011 

0.042 

0.162 

40 

10 

4.0 

0.004 

0.009 

0.980 



4.5 

0.001 

0.002 

0.440 



5.0 

0.000 

0.002 

0.034 


Table 5.7 Temporal moments for different values of k 


k 

Analytical Moments 

Computed Moments 

Zeroth 

First 

Second 

Third 

Zeroth 

First 

Second 

Third 

0.03 

1.000 

6.436 

14.852 

70.880 

1.000 

6.436 

14.852 

70.880 

0.3 

1.000 

6.436 

2.735 

2.365 

1.000 

6.436 

2.735 

2.365 

3.0 

1.000 

6.436 

1.523 

0.965 

1.000 

6.436 

1.524 

0.964 
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Table 5.8 Temporal moments with various values of a 


a 

Analytical Moments 


Computed Moments 

Zeroth 

First 

Second 

Third 

Zeroth 

First 

Second 

Third 

“M04 

1.000 

6.436 

5.287 

24.774 

1.000 

6.436 

5.287 

24.774 


1.000 

6.436 

2.735 

2.365 

1.000 

6.436 

2.635 

2.635 

0.4 

1.000 

6.436 

2.480 

1.871 

1.000 

6.436 

2.480 

1.871 , 
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Figure 5.3 Concentration profile at quarter length of the 

column 



Figure 5.4 Concentration profile at half length of the 

column 
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tration 








Figure 5.9 Effect of k on skewness of the break-through 

curve 



Figure 5.10 Effect a of on the skewness of the break- 
through curve 
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Figure 5.13 Concentration profile at mean arrival time 
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Figure 5.14 Comparision with field observations 
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